IUCAA-22/2001, KSUPT-01/3 May 2001 



Window Function for Non-Circular Beam CMB Anisotropy Experiment 

Tarun Souradeep 1,2 and Bharat Ratra 1 



ABSTRACT 

We develop computationally rapid methods to compute the window function for a 
cosmic microwave background anisotropy experiment with a non-circular beam which 
scans over large angles on the sky. To concretely illustrate these methods we compute 
the window function for the Python V experiment which scans over large angles on the 
sky with an elliptical Gaussian beam. 

Subject headings: cosmic microwave background — cosmology: theory — methods: analy- 
tical — methods: data analysis 



1. Introduction 

Cosmic microwave background (CMB) anisotropy measurements are becoming an increasingly 
powerful tool for testing cosmogonies and constraining cosmological parameters. See, e.g., Subrah- 
manyan et al. (2000), Romeo et al. (2001), Dawson et al. (2000), and Padin et al. (2001) for 
recent CMB anisotropy observations, and, e.g., Ratra et al. (1997), Gorski et al. (1998), Rocha et 
al. (1999), Gawiser & Silk (2000), Knox & Page (2000), Douspis et al. (2001), and Podariu et al. 
(2001) for discussions of constraints on models from the CMB anisotropy data. 

Conventionally, the CMB temperature, T(-y), is expressed as a function of angular position, 
7 = (6, (p), on the sky via the spherical harmonic decomposition, 



T ^ =EE a emY em (l) • (1) 



£=0 m=-l 



The CMB spatial anisotropy in a Gaussian model 3 is completely specified by its angular two-point 
correlation function C(-y,-y'), = (T(~f)T(~f')), between directions 7 and 7' on the sky. In most 
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3 The simplest inflation models predict a Gaussian CMB anisotropy (see, e.g., Fischler, Ratra, & Susskind 1985) 
on all but the smallest angular scales. CMB anisotropy observations on quarter degree and larger angular scales 
appear to be Gaussian (see, e.g., Mukherjee, Hobson, & Lasenby 2000; Aghanim, Forni, & Bouchet 2001; Phillips & 
Kogut 2001; Park et al. 2001; Wu et al. 2001b, also see Podariu et al. 2001). 
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theoretical models the predicted fluctuations are statistically isotropic, (7(7,7') = (7(7-7'). The 
fluctuations can then be characterized solely by the angular spectrum Cg, defined in terms of the 
ensemble average, 

(ae m ae>m'*} = C(8u>S mjn i , (2) 
and related to the correlation function through 

<r(7) r(V)> = f; Q P^(7 • V) , (3) 

where is a Legendre polynomial. 

Typically, a CMB anisotropy experiment probes a range of angular scales characterized by 
a window function W^-y, 7'). 4 To utilize the full information in the data one must use model 
anisotropy spectra (Q's) defined over this range of angular scales. Such theoretical spectra are 
parameterized by cosmological parameters such as f-o, h, and in these models 5 , and by the 
spectrum of quantum fluctuations generated during inflation. 

To use model GYs in conjunction with CMB anisotropy data to estimate cosmological parame- 
ters one must be able to carefully model the CMB anisotropy experiment, i.e., accurately compute 
the window function We. Given such a model of an experiment and a family of Ces, one may opti- 
mize the fit to the data from the experiment either by using an (approximate) x 2 technique (e.g., 
Ganga, Ratra, & Sugiyama 1996; Bond, Jaffe, & Knox 2000; Knox 1999; Rocha 1999; Lineweaver 
2001; Dodelson 2000; Tegmark Sz Zaldarriaga 2000) or by using an exact maximum likelihood 
technique (e.g., Gorski et al. 1995; Ganga et al. 1997, 1998; Ratra et al. 1998, 1999). 

Current CMB anisotropy data is of significantly higher quality than data that was available 
just a few years ago. As a consequence an accurate model of an experimental We must now account 
for effects that were ignored for earlier experiments. In this paper we develop computationally rapid 
methods that account for the non-circularity of the beam in a CMB anisotropy experiment window 
function at large angular separations where the curvature of the sky cannot be ignored. This must 
be accounted for in an experiment like Python V (Coble et al. 1999) which has an elliptical beam 
and samples a large enough area of the sky to prejudice use of the flat-sky approximation. While 
we focus here on Python V as a concrete illustrative example, our techniques are easily generalized 
to more complex cases (e.g., arbitrary beam shape, beam rotation, non-Gaussianity of the beam, 
etc.). Wu et al. (2001a) develop an alternate method to deal with the asymmetric beam of the 
MAXIMA- 1 experiment. 



4 See www.phys.ksu.edu/~tarun/CMBwindows/wincomb/wincomb_tf.html for a discussion and tabulation of zero- 
lag window functions. 

5 Here f_o is the nonrelativistic-mass density parameter, h is the Hubble parameter in units of 100 km s _1 Mpc -1 , 
and fis is the baryonic-mass density parameter. 
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In §2 we describe the general formalism for computing the window function. In the §§3- 
5 we develop specific computationally rapid methods for computing the window function in three 
different cases. The flat-sky approximation window function computation is covered in §3. In §4 we 
develop a general method, based on Wigner rotation functions, for computing the window function 
on a sphere, and describe how to numerically implement this scheme. In §5 we specialize to the case 
of an experiment like Python V where long scans are performed at constant elevation, and provide 
a computationally rapid method for evaluating the exact window function. Approximate window 
functions obtained using the flat-sky approximation and from retaining only the first few terms in a 
perturbation expansion (in non-circularity about a circular beam) of the Wigner rotation functions 
method (hereafter Wigner method) are compared to the exact window function in §6. We conclude 
in §7. In appendix A we describe our parameterization and normalization of an elliptical Gaussian 
beam, and also record analytic expressions for its Fourier and spherical harmonic transforms. 



2. Window Function Computation Formalism 

Owing to the finite angular resolution of a CMB anisotropy experiment, the temperature 
"measured" by the experiment at point 7i on the sky is 

f( 7i ) = j ^ 7 B( 7i ,7)r(7). (4) 

Here -B( 7i , 7 ) is the beam function that characterizes the angular dependence of the sensitivity of 
the apparatus around the pointing direction 7i . 

CMB anisotropy experiments that use a differencing or modulation scheme measure the dif- 
ference in temperature between different points on the sky. The measured CMB temperature 
anisotropy in any differencing scheme (labelled below by index n) can be expressed as a weighted 
linear combination, 

A (n) ( 7 J = J ^ 7 ^ (n) (7)T(7), (5) 

where w\ a \-y) are the weight functions. In our discussion of Python V below n corresponds to 
the harmonic number of cosine modulation weight functions. Instead of weights ^"^(7) that are 
continuous functions, for Python V the weights are discrete and characterize the sensitivity 
at 7jj, the j th point on the discretized chopper cycle around the pointing direction 7i . In this case 
the integral in eq. (5) is replaced by a summation and we have 

A^(7j = X;^; ) r(7 ij ) (6) 
3=1 

where N c is the number of points in the discretized chopper cycle. 

The complete window function for modulation pair (n,m), W^ n,m ^(7j, 7j), accounts for the 
effects of both the beam function and the differencing or modulation scheme of the experiment and 
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is defined through the theoretical model covariance matrix 

C (n,m) = <A (n) (7 . )[A (m) (7 . )r> = £ (^tll Q ^( 7j , 7 .) . (7) 

e=o 

It proves convenient to distinguish between the dependence of the window function on the finite 
angular resolution of the experimental apparatus (the beam function) and the dependence of the 
window function on the differencing scheme adopted for the experiment. Using eq. (6), the complete 
theoretical covariance matrix element C^' m ^ between pixels i and j on the sky can be expressed 
as a weighted linear sum of single beam correlation functions via 

N c N c 

ctr ] -EE -l n) <r(7*)r(7 iff )> • (8) 

p=l q=l 

We call the single beam correlation function (7^(7,7'), = (T( j j)T( j j')) , an elementary correlation 
function. The elementary correlation function does not depend on the differencing scheme used in 
the experiment but does depend on the beam function. We use to define what we call the 
(single beam) elementary window function via 

^ (e) (7„7 J ) = E^^Q< ) (7 i ,7 J ). 0) 
e=o 

(e) 

Wp depends on the beam function of the experiment but not on the differencing strategy used. 
Using eqs. (8) and (9) the complete window function may be expressed as a weighted linear sum 
of elementary window functions via 

<' m) (7„ 7,) = ee-s;M; } < 7*) • (io) 

p=l q=l 

Using eqs. (3), (4), and (9) the elementary window function may be expressed as 
Wt\l^ 7j ) = J ^ 7 dn r B( 7t , 7 )i?(7„ 7)^(7 ■ V) 



4tt 



2£ + 



T E beMlbemdj)}*, (11) 



where 



hmili) = J dfi 7 B (7i, 7)^(7) (12) 
is the spherical harmonic transform of the beam function pointing at 7i . 

For some experiments the beam function is accurately circularly symmetric about the pointing 
direction, i.e., £(7^,7) = -6(7^ -7). This allows a great simplification 6 since the beam function can 



6 When the beam is pointing towards a Pole, 7 P , the coefficients in the spherical harmonic expansion of a circular 
beam function are bi m (-y P ) = Bt\J (2£ + l)/4ir 8 m o- Here 5 m o is a Kronecker delta function and implies that the only 
non-zero spherical harmonic coefficients are those that result in a circularly-symmetric beam function. 
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then be represented as B{^f i -7) = (4ir) 1 (2^ + 1) Bg Pi (7, -7). Consequently, for a circularly 

symmetric beam function it is straightforward to derive the usual expression 



In addition to the single evaluation of the Legendre transform of the beam, Bg, computation of the 
elementary window function for an experiment with a circular beam simply involves computing Pi 
for all I < £ max using the stable upward recursion relation for each distinct pixel pair separation. 
This is computationally inexpensive (at least by a factor ~ £ max ) compared to the computation of 
the elementary window function for an experiment with arbitrary beam shape (eq. [11]). 

In the next three sections we discuss three cases where one may compute the window function 
for an experiment with a non-circular beam in a time comparable to or just a factor of few larger 
than that for the same experiment assuming a circular beam. First, we consider the flat-sky 
approximation which is accurate if the experiment has a compact beam and the pixels are not 
separated by large angles (more precisely, the separation must be significantly less than a radian). 
This has been used previously for a number of experiments, including Python V (Coble et al. 1999; 
Coble 1999), MSAM (Coble 1999), and MAXIMA-1 (Wu et al. 2001a). Next, we develop a very 
general Wigner method that fully accounts for the curvature of the sky and accounts for the non- 
circularity of the beam in a perturbative expansion about a circular beam. In a case like Python 
V where the non-circularity of the beam is not large, the first few terms in the non-circularity 
perturbation expansion provide sufficient accuracy. In this case the Wigner method allows one to 
compute the window function in a time a factor of a few larger than that for the corresponding 
circular beam case. Finally, for an experiment that scans at constant elevation (such as Python 
V), with pixels lying on a small number of elevations, it is possible to implement a slightly different 
Wigner method that allows rapid computation of the exact window function for an arbitrary beam 
shape. 



If an experiment scans a small enough patch of the sky it is computationally advantageous 
to work in the flat tangent plane (rather than on the sphere) and make use of Fourier transforms 
(rather than spherical harmonic transforms) when modeling the experiment. See, for instance, 
Bond & Efstathiou (1987) and Coble (1999) for discussions. We may then transform from 7 to 
coordinates in a locally flat patch, u, and use a two-dimensional Fourier transform approximation 
to the spherical harmonic transform. For instance, 



W r i e) (7i,7i) = 5|P/( 7i -7 J -)- 



(13) 



3. Window Function in the Flat-Sky Approximation 




(14) 



Here u, = (to, </>), are polar coordinates in the neighborhood of the North Pole 7 P in the patch on 
the sky, i.e., 7 = 7 p + cj, or u> = (u>i, oj 2 ) = (ujcoscf), uismcj)), where uj = 2 sin(0/2) = I7 — 7p|, < 
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lo < 2, where 9 is the colatitude (e.g., Bond & Efstathiou 1987). In the small-angle approximation, 
it may be shown that the ensemble average of the Fourier transform of the temperature 

(T(k)T*(k')) = (2^) 2 C fc <5( 2 )(k - k') , (15) 

where k = |k| and Ck is the CMB anisotropy power spectrum. This flat-sky analog of eq. (2) is 
obtained assuming statistical homogeneity on the flat sky. The angular correlation function can 
then be expressed as 

(T(u l )T(u j )) = ±- J^dkkMklui-UjDCk, (16) 

where Jo is the zeroth order Bessel function of the first kind. Comparing this expression to that 
in eq. (3) in the small angular separation and large £ regime where Pi(cos9) — > Jo([£ + l/2]6), we 
arrive at the correspondence k ~ £ + 1/2 between the radial wavenumber on the flat sky k and the 
spherical multipole £. 

For an experiment for which the flat-sky approximation is valid the "measured" temperature 
(see eq. [4]) is 

where -B(k) is the Fourier transform of the beam function pointing at the origin u = of the 
local flat coordinate patch. The rotation operator Rj (which rotates k by an angle q{) accounts 
for a possible rotation of the telescope beam pointing at u>j relative to the telescope beam pointing 
at the origin. In addition to the case when the telescope physically rotates around its axis as it 
moves from one pointing direction to another, non-zero values of Qi can arise when the telescope is 
not located at a Pole and also when a single flat-sky coordinate system is set up on a patch large 
enough for sky curvature to be important. In the latter two cases this rotation is important only 
in the regime where one expects the flat-sky approximation to be poor. 

The modulated temperature A( n )(7j) expressed in terms of T(u)ij) is given by eq. (6). The 
window function wj> n ' m ^ (7, , 7^ ) is defined through the covariance matrix by 

c (n, m) = £ (2£+l) Ce W (n,m) { ^ ^ cT^kC, W^\ lh 7 ,) , (18) 

where we identify the flat space radial wavenumber k with £+ 1/2. The complete window function 
W^ n ' m ^(7j, -fj) may be expressed in terms of elementary window functions wjf^ by using eq. (10). 

Using eqs. (9) and (17), the expression for the elementary window function in the flat-sky 
approximation is 

Wi e \u h Uj ) = fj ^ e^^-^)S(R t [k])i?*(R, [k]) . (19) 
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Further analytical manipulations are needed to derive an expression suitable for numerical eval- 
uation. Without loss of generality, we transform to a new flat coordinate system with u>j as the 
origin: uj\ = (0,0) and u'j = {lo'^uj 1 ^). In what follows we drop the prime on the new coordinates. 
For small angular separations, to = — Uj\ <C 1, we have 

W^(u t , Uj ) = r ^ cos [ku;co S (<p k -a)} B(k) B*(Rj[k]) , (20) 
Jo Z7T 

where we have defined a = tan~ x [ui2 / Vi) • 

For an elliptical Gaussian beam function (eq. [Al]), an analytic expression for B(k) is given 
in eq. (A7). Using this in eq. (20) (and ignoring rotations, Qj = 0) we obtain 

W { °\u h Uj) = I — cos[kucos{(f> k - a)] exp [-fcVi 2 (l + esm 2 (j) k )] , (21) 
Jo 71 

where u\ and o"2 are the beamwidths (in radians) along the major and minor axis of the elliptical 
Gaussian beam. It is useful to define the RMS beamwidth ctrms> = y (o'l 2 + °2 2 ) /2, and a measure 
of the non-circularity of the elliptical beam e, = (02/cxi) 2 — 1- 



The expression for the flat-sky window function for an elliptical Gaussian beam given in eq. (21) 
can be readily evaluated numerically. It can also be expressed analytically as an infinite series 
expansion, wjf^ = ^2'^L (— l) ra ^wjf \ in powers of an anisotropy or non-circularity parameter 
b k = e(aik) 2 = k 2 {a\ — af). The series is perturbative (b k < 1) up to multipoles t"z\(j\ — cr 2 | -1 / 2 , 
which can be considerably larger than the inverse of the RMS beamwidth <trms- As expected, the 
first (n = 0) term in this series corresponds to the circular beam function result (a\ = 02 = c), 

2 



For a non-circular Gaussian beam function the next few terms are 



(22) 



J\{kuj) 



kuj 
J 2 (kuJ 



— J2(kuj) sin 2 a 
Jz(kuj 



6- 



sin a + J/±(koj) sin a 



2 [" (koj) 2 " kuj 

where J m (x) is the m th order Bessel function of the first kind. At arbitrary order n the term is 



(23) 
(24) 



¥r(n + l/2) _ k 2 ai 2 ^ (sinQ) 2m (cosa) 2 ( n - m ) 



r(n+l) 



E 

m=0 



r(n — m + 1) r(m + 1) 



x jF 2 



- + m, -, 1 + n 



(A; a; 



,21 



(25) 



where T is the Euler gamma function and 1F2 is a generalized hypergeometric function. Figure 1 
shows contour plots of the zeroth order (eq. [22] ) and first order (eq. [23] ) terms in the non-circularity 
expansion of the flat-sky window function. These are computed for parameter values characterizing 
the Python V experiment (see §6 below for details of the experiment). 
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4. Wigner Method Window Function 



If an experiment takes data over a large enough area of the sky the formalism developed in 
the previous section, based on the flat-sky approximation, cannot be used to compute the window 
function. In this section we develop a general method for computing the window function for an 
experiment with an arbitrary beam shape that collects data from a large area on the sky. 

For pointing direction = (9i, fa), and vector in the beam 7, = (9, fa), the beam function 
may be expanded in a spherical harmonic decomposition, 



Here the expansion coefficients, bi m , are given by eq. (12). 

For ease of computation it is convenient to rotate to a new coordinate system in which the new 
X3 axis lies along the pointing direction y i . This is accomplished by first rotating the coordinate 
system around the x 3 axis by fa and then rotating around the new x' 2 axis by 0j. Then 



where the last step follows from the fact that B is a scalar. 

If the experiment is not located at a Pole, the beam of a telescope which does not rotate 
around its beam axis will, nevertheless, appear to "rotate" around the beam axis with respect to 
the local azimuth and declination directions. We hence allow for a rotation gi of the beam around 
the beam axis relative to a parallel transport of the beam on the sky from the Pole 7p, = X3, to 
the pointing direction 7^. The Python V experiment was located at the South Pole and hence has 
Qi = 0. The rotation Qi can also account for non-circular beam function cases where the telescope 
rotates around its axis as it moves from one pointing direction to another (e.g., one mounted on a 
satellite) . 

The rotations of the previous two paragraphs correspond to Euler angles a. = fa, (3 = 9i, and 
7 = Qi in the notation of Scheme A of §1.4.1 of Varshalovich, Moskalev, & Khersonskii (1988, 
hereafter VMK). From eq. (1) of their §5.5.1 we have 




(26) 



£=0 m=-i 




(27) 



£'=0 m'=-l> 



11 




(28) 




is a Wigner D-function corresponding to the rotation of the beam. From eqs. (28) 



and (27) we have 



I 




(29) 



rn ' 



■> = -l 
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We focus on the elementary window function, i.e., ignore the modulation and consider the 
window function for two points f y i and "y ■ (see eq. [11]). Using the usual decomposition for the 
Legendre polynomial in terms of spherical harmonics, the fact that B is real, and eqs. (28) and 
(29), we find 



e 



' -—l m"=-e 



E [ D lm'(^ Oi, Si)}* Dim»(h, Oj, Qj)- (30) 



X 

m=-i 



When the beam function is circularly symmetric 6^ m (xg) = 5 m oB£ y/2£ + 1. Using _D^ O (0, <fi, g) = 
Y£ m (6, (ft), it is straightforward to establish that in this case eq. (30) reduces to the usual expression 
given in eq. (13). 

Using the addition theorem for Wigner /^-functions (eqs. [2] of §4.4, [1] of §4.3, and [5] and 
[6] of §4.7.2 of VMK), we reduce eq. (30) to the simpler form, 



47T 



I 



WT(7i, -Tj) = 2eTl ? ? [W(x^)]*W(x^ m »(«-ft,7,/3 + ^), (31) 



m'=—t m"=—l 



where 



COS7 = cos^j cosdj + sin^j sin0j cos(0j — 4>j) = 7j • 7^ 

cota = — cos#j cot((f)i — + sin6>; cot^j csc(4>i — <pj) (32) 
cot/3 = —cosOj cot(4>i — (fij) + cotOi sin9j csc(<^ — <pj). 



For large values of £ it is computationally expensive to evaluate the entire m! and m" sum 
in eq. (31). However, for a smooth, mildly non-circular beam function (defined precisely below), 
restricting the summation to a few low values of m' and m" results in a good approximation. A 
smooth (i.e., not sharply peaked) beam function results in falling off with increasing I. At 
large I and small 7, D l mm ,{a, 7, f5) e~ ima ~ im ' I 3 J\ rn _ ml \{[£ + l/2] 7 ). Thus the D e m , m „ term in 
eq. (31) strongly suppresses, as (l^ 111 ~ m the contribution from off-diagonal m! 7^ m" terms. In 
addition, mild non-circularity requires that at each value of £ the ratio \b( m /beo\ decrease rapidly 
with increasing \m\. Hence the products b^ m ,b( m ", as functions of £, are ordered in magnitude 
and fall off as one goes to higher values of \m'\ + \m"\. For Python V where the deviation from 
circularity is small, retaining the first non-zero order term is sufficient for computing an accurate 
covariance function. 7 



7 By a sufficiently accurate covariance function we mean that the maximum likelihood analysis results are not 
significantly affected by higher order corrections. This criterion of sufficient accuracy therefore depends on the level 
of noise in the experiment, which is encoded in the noise covariance matrix of the experiment. 
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We now derive explicit expressions for the first few leading order terms in eq. (31) for the 
specific case of an elliptical beam function. For an elliptical beam function, symmetry dictates that 
be m (x' 3 ) = for odd m (see appendix A). In what follows beam rotations are set to zero (gi = 0), 
but it is straightforward to restore them from the complete expression given above. The zeroth 
order term contains Dq , the first order term has four contributors (D^ 2 , Dq _ 2 , D 20 , and D £ _ 20 ), 
and the second order term has eight contributors (D 22 , ^2 -n ^-2 2> -^-2 -2' ^04' -4> -^40> anc ^ 
/>' ,oi- " ' " 

After a significant amount of algebraic manipulations (that use equations from §§4.3, 4.4, 4.8, 
and 4.17 of VMK, familiar properties of Legendre polynomials and modified Bessel functions, and 
the reality condition on the beam, eq. [A8]), we find a series expansion of the elementary window 
function, 



21 + 1 



2 4(7) + 2M4)M4) {cos(2q) + cos(2/3)} <&( 7 ) 
+2 [b e2 (x' 3 )] 2 [cos(2a + 2/3)42(7) + (-l)^cos(2a - 2/?)4(tt - 7) 



+26^o(x / 3 )6 £ 4(x / 3) {cos(4a) + cos(4/?)} ^ 4 ( 7 ) + • • • 



(33) 



Here the angles a, /3, and 7 are defined in eqs. (32), and the 4m' ' s are the usual Wigner d- 
functions of angular momentum theory (e.g., §4.3 of VMK) related to the Wigner D-functions 
through D e mm ,(a,"{,f3) = e ~ tma d e mm , (7) e ~ tm ' ^ . More precisely, 



4(7) 
4(7) 

4(7) 



P^(cos7), 



+ 1) 



(£-!){£ + 2) 



Pj?(cos7) + 



2cos7 



1 



(*-!)(* + 2) 



-4|^^)+^ + l) 



1 + COS7 



P^(cos7) 



P £ '(cos7), 



4(1 — COS7) 



(£-!)£(£ + !)(* + 2) 



C0S7 



1 + C0S7 



+ €(£ + 11 



P/(cos7), 



(34) 



(-1)* 



(l-W + 2) 



_ 4( '2 + cos7^ + 



C0S7 



P^(cos7) 



+ 



4(7) = -^+1 



4(-l)^ +1 (l + COS7) I" /2 + COS7 
(£- 1)^(£ + l)(€ + 2) [~ V1-COS7 

6(1 + 3cos 2 7) 



+ £{£ + !) 



P/(cos7), 



'(£-4)1 
(£ + 4)1 



sin 2 7 



- £(£ + 1) 



Pi(cos-y) 



+8. 



4)! 



(£ + 4) 



C0S7 



3(1 + cos 2 7) 



sin 2 7 



- £(£ + 1) 



P/(cos7) , 



where P^' = dPe(x)/dx. Recursion relations in the indices £, m and m' (see VMK) can be used to 
compute d\ 1 , for larger values of m and m! . 
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Evaluating the first few terms of the Wigner method expansion in eq. (33) involves computing 
Pi and PL The first derivative P[ can be readily computed in terms of Pi and Pi-\ during the 
generation of Pi using the upward recursion relation. Hence the computational cost of evaluating 
for a mildly non-circular beam function (using eq. [33]) is only a factor of a few larger than 
that for a circular beam function. 

Figure 2 shows, as a function of ^ormS; the six leading order \bi m (x' 3 ) 6£ m '(xg)| coefficients 
(with m,m' = 0,2,4) of the expansion of eq. (31) (see eq. [33]) for the Python V experiment. 
These are also the leading order contributors to the zero-lag elementary window function. Note 
that the curves do not cross at any £, i.e., the ordering of the coefficients is maintained for all £, and 
at large £ (past the peak) higher order terms fall off more rapidly with i. Hence this perturbation 
expansion is an efficient scheme for computing non-circularity corrections. The non-circularity 
corrections peak at angular scales smaller than the RMS beamwidth orms- Thus non-circularity 
corrections are not that important for an elementary window function (a single beam experiment) 
but can have a significant effect on the complete window function for a modulated experiment if 
the modulation scheme results in sensitivity to the ^orms > 1 regime. 

Figure 3 shows contour plots of the isotropic term c?q (7) and the leading order correction term 
{cos(2a) + cos(2/3)} d l m {l) in the perturbation expansion of the Wigner method elementary window 
function (see eq. [33]) for the Python V experiment. These are plotted for multipole £ = 100 chosen 
so that ^orms ~ 1) which is where the non-circularity correction starts becoming significant. Close 
to the center of the plots this non-circularity correction is larger along the major and minor axes 
of the elliptical beam. The correction term falls off with increasing pixel separation, suggesting 
that the circular beam function approximation is good for sufficiently large separations. Also, for 
a modestly non-circular beam function this implies that higher order terms in the Wigner method 
perturbation expansion need be retained only for close pixel pairs and one can hence truncate the 
summation in eq. (31) at lower orders for progressively more widely separated pixel pairs. 



5. Constant Elevation Scan Window Function 

For an experiment, such as Python V, that scans at constant elevation 8 , it is possible to derive 
another expression for the window function which does not require use of the approximation of the 
previous section (the truncation of the m' and m" series). 

We follow the initial development of the previous section and transform to a new coordinate 
system by rotating around the X3 axis by fa, where j { = (0{, fa) is the pointing direction. The 
Euler angles of this rotation are a = fa, (3 = 0, and 7 = Qi. As described in the previous section, 
Qi represents a relative rotation of the beam around its axis. For an experiment located at one of 



Constant elevation here refers most generally to any set of parallel circles on the sky. 
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the Poles ft = 0. In general, we have 

t 

bimhi) = E ^'(^s)^^.. °> (35) 

m'=-l 

From §4.16 of VMK we find £>^ m ,(<fc 0, ft) = e -™(*i+«) s mm ,. Equation (35) then implies 
be m (9i, <f>i) = e- im ^ i+ ^b em (6i, 0) and so eq. (30) reduces to 

Wt\li, lj) = 2^ E e-M^J+te-ft)] 6| m (^, 0) 6, m (%, 0). (36) 

m=— i 

Window functions between pixels lying on a few constant elevation lines can be rapidly computed 
by using eq. (36) and pre-computed bi m (9i, 0)'s. 

For an experiment like Python V whose beam function has the symmetry B(9i, 0; 6, 4>) = 
B(6i, 0; 9, —(/)), it may be shown that eqs. (12) and (26) and the fact that the beam function B is 
real implies that be m (9i, 0) is real. In this case eq. (36) may be re-expressed as 

Wt\li, 7,) = ^77 E cos " + " be m (9i, 0) b tm {9j, 0). (37) 

m=—£ 

Here bi m {9, 0) is defined in the usual way through eq. (12). For the Python V experiment, which 
chops at constant elevation and where the pixels lie on a relatively small number of elevations 9 , it 
is possible to pre-compute and store the be m (9, 0) at all elevations 9. Given these pre-computed 
bimS, the elementary window function can be computed very rapidly by using eq. (37). 



6. Comparison of Approximate and Exact Python V Covariance Matrices 

Python V is a CMB anisotropy experiment which performs wide-angle scans with a non-circular 
beam. In this section we compare Python V window functions and theoretical covariance matrices 
computed in the flat-sky approximation and in the Wigner method perturbative approximation as 
well as in the exact method. 

Python V observations were made at 37-45 GHz. Two regions of the sky were observed: the 
main field, a rectangle 7?5 in declination (<5 = —52° to — 45?4) by 67?7 in azimuth (centered on 
a = 23^18); and another rectangular patch 3° in declination (5 = —63° to —60°) by 30° in azimuth 
(centered on a = 3^0). For detailed descriptions of the experiment and data see Coble et al. (1999) 
and Coble (1999). 

The Python V beam function is well described by an elliptical Gaussian with FWHM beam- 
widths of l?02l°;°f in elevation and 0?91±g|g? in azimuth (one standard deviation uncertainties). 



'The 690 sky pixels of Python V lie on 11 distinct elevations. 
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The Python V beam function is a compact elliptical Gaussian, eq. (Al), with <j\ = 0.0076 and 
(7 2 = 0.0067 as the nominal beamwidths in radians. Python V uses a constant-elevation smooth scan 
sampling strategy around every pixel on the sky, with a chopper throw (end to end) of $ c = 17?06. 
Each chopper cycle consists of 128 time samples suitably modulated in time to correspond to the 
spatial modulations described below. To compute the window function we only need to know 
the final spatial modulation strategy adopted. See Coble (1999) for a more detailed discussion of 
these procedures. The constant elevation scans are discretely resampled in space (with nine-fold 
oversampling) at N c = 567 equi-spaced points (6i, labelled by integers p = 1,2, . . . ,N C along the 
chopper cycle around the pixel 'y i = (f>i). The azimuth (f>? = 4>i + (p— l)A<f> + (<J>„ — 3> c /2), where 
A(j) = 3> c / (N c — 1) is the spacing between points and 3> = 0?58 accounts for the offset between 
the azimuth of the pixel, fa, and the center of scan. The scans are modulated using the first eight 
cosine harmonics of the chopper cycle (hereafter modulations 1 to 8). All modulations, other than 
the first, are apodized by a Hann window to reduce ringing in multipole space and down weight 
data taken during chopper turnaround. For the modulated Python V scans the weight functions 
(see eq. [6]) are 10 

< ^ 2M (m) 
( m ) _ IMp 

where 

MW = cos(2ttZ p ) 

f_1 \m+l 

M W = - — '- cos(27rmZ p ) [1 - cos(2ttZ p )} , (m > 1) (39) 

with Z p = (p — 1)/(N C — 1). These weights are used to obtain the Python V complete window 
function W^ n ' m \ , y i , ~fj) (between modulation m of the scan around sky pixel 7j and modulation 
n of the scan around sky pixel *y J -) in terms of the elementary window functions (see eq. [10]). 

The Mp m ^ are identical for all pixels hence is independent of pixel index i. With identical 
constant-elevation chops around every pixel, eq. (10) can be re-expressed in a computationally more 
efficient form. In this case the elementary window function wj^ depends only on (fl? — 4> q - (at fixed 

6i and 6j). This reduces the number of separations at which is needed and thus speeds up 

the computation. In this case we may re-express the window function as 



N c -1 N c -qi 

^ (n ' m) (7„7,o= y, y. 



91 = 1 92 = 1 



10 The following expressions correct a typographical error in the corresponding expressions in Coble et al. (1999) 
and Coble (1999). 
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N c -1 
92=1 

Here the three arguments of the elementary window functions (the three between two direc- 
tions) are the colatitudes of the two directions and the difference in their azimuth angles. Pre- 
computing the second summations over the weight products, vop 1 ^ = Y^=\ P w i?q+p) w jq^ ^ or tnc 
set of all the modulation pairs speeds up the evaluation of the right hand side of eq. (40) . 

Figure 4 shows the eight equal-modulation exact Python V zero-lag complete window functions 
at the two extreme values of the declination, S = —63° and 5 = — 45?4. Also plotted for comparison 
is the zero-lag elementary window function at 5 = —63°. The effect of beam function non-circularity 
is more pronounced for the modulated (complete) window functions since these have peak sensitivity 
at multipoles well beyond the inverse beamwidth, which is where the non-circularity corrections 
start to become important (see Fig. 2). 

While it is of interest to estimate the accuracy of window functions computed in various 
approximations, the accuracy of computed covariance matrix elements is of much greater relevance 
since these directly determine the accuracy of cosmological results extracted from CMB anisotropy 
data. This significantly extends the range of usefulness of approximate window functions. First, the 
sum over I in the definition of the theoretical model covariance matrix, eq. (7), averages over and 
hence reduces the significance of deviations between the approximate and exact window functions. 
For example, an approximate window function that has large deviations from the exact window 
function only in regimes where they oscillate in i may still result in an accurate covariance matrix. 
Second, errors in the window function are unimportant when the corresponding covariance matrix 
elements are small (subdominant). For example, while the flat-sky approximation window function 
is very inaccurate for a widely separated pixel pair, the covariance matrix element for such a pixel 
pair is subdominant and thus cannot significantly influence the results from a maximum likelihood 
analysis of the data. Another important consideration is the level of noise in the experiment. 
Since the inverse of the sum of the theoretical model and noise covariance matrices, (Ct + Cat) -1 , 
determines the results of the maximum likelihood analysis, the model covariance matrix needs to 
be computed to higher accuracy for an experiment with lower noise. 

Equation (7) defines the model covariance matrix element C^™ - '™^ in terms of the window func- 
tion and the model CMB anisotropy power spectrum. To compare covariance matrices computed 
using different approximations therefore requires choice of a model Cg. We use the flat bandpower 
spectrum, Cg oc l/£(£ + 1), for this purpose in what follows. With this choice of power spectrum, 
the quantity 

*f' m) (7i,7,-) = £ 7^)<' m) (7. 7,) , (41) 

is a measure of the cumulative build up of the corresponding covariance matrix element as one 
progresses with the sum over multipole t in eq. (7). Up to a multiplicative normalization constant, 
-^ n '"0(7i>7j) conver g es to the covariance matrix C^' m ^ as t — ► oo. Unlike Wg, Fg measures ^-space 
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information that may be directly used to estimate the accuracy of the computed covariance matrix. 
It also allows one to determine the value of t to which one must compute to achieve a desired 
accuracy. 

Figure 5 compares Python V flat bandpower covariance matrix elements computed using the 
different approximations developed above. The upper panels, a) and b), show comparisons of zero- 
lag covariance matrix elements. 11 At zero-lag the flat-sky approximation is more accurate than 
the zeroth order (circular beam) and first order (wig20) Wigner method perturbation expansion 
approximations. The Wigner method gets progressively more efficient (i.e., one needs to retain 
fewer terms in the perturbation expansion series to achieve the desired accuracy) for non-zero-lag 
covariance matrix elements between increasingly separated pixels. For covariance matrix elements 
between very widely separated pixels even the circular beam (zeroth order) approximation suffices. 
The Wigner method is more accurate for lower modulations which probe lower values of I where 
beam non-circularity corrections are smaller (see Fig. 5). The flat-sky approximation is accurate at 
small separations (e.g., at zero-lag). For Python V the flat-sky approximation is more accurate for 
pixels separated in azimuth than for pixels separated in declination. The flat-sky approximation 
works better for higher modulations which probe larger values of I. The F( curves for the circular 
beam approximation (wigOO) in the lower left panel c) also highlights the possible pitfall of not 
computing to large enough £. Here F^j does not converge to Cjj until £ > 500 ~ 4cr^^ Is . 



7. Conclusion 

We develop computationally rapid methods to compute the window function for a long-scan 
arbitrary beam shape CMB anisotropy experiment. We use these methods to compute the window 
function for the elliptical Gaussian beam Python V experiment. 

It proves convenient to separate effects due to the modulation scheme adopted and the shape 
of the beam function by expressing the complete window function as a weighted sum of single-beam 
elementary window functions, eq. (10). 

Using eq. (11) to obtain exact elementary window function for a non-circular beam experiment 
requires accurate computation of the spherical harmonic transform of the beam function at each 
pointing direction. For an experiment with a large number of pointing directions this is computa- 
tionally prohibitive. For instance, Python V has 690 pixels and the scan around each is resampled 
at 567 points which results in ~ 0.4 million pointing directions. Fortunately, the 690 pixels lie on 
only 11 distinct elevations and the scans are performed at constant elevation. Hence pre-computing 
and storing the spherical harmonic beam function transforms at these 11 declinations allows for 
rapid computation of the exact Python V window function. 



11 Modulated zero-lag complete window functions receive contributions from non-zero-lag elementary window func- 
tions between pixels separated by as much as the chopper throw 3> c = 17? 06. 
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In the absence of such a simplification due to a symmetry, the Wigner method perturbation 
expansion scheme allows an accurate computation of the window function with computational ef- 
fort within a factor of a few of the corresponding computation for the case of a circular beam 
function. This factor depends on the order to which the perturbative expansion (about the circular 
beam approximation) needs to be developed to achieve the desired accuracy. In this implementa- 
tion, the Wigner method requires pre-computation and storage of one spherical harmonic trans- 
form of the beam pointing at a Pole, be m (~fp). If the beam is mildly non-circular then for all £ 
\be m ('Jp)\/\beo(lp)\ fehs off rapidly with increasing \m\, allowing for a fast and accurate transform 
and simpler storage. In the Appendix we record a semi-analytic expression for the beam function 
transform for a compact elliptical Gaussian beam. In mildly non-circular cases, such as Python 
V, |^£m(7p)|/|^£o(7p)l rapidly with increasing \m\ and the first order Wigner method is suffi- 
ciently accurate. We also develop a flat-sky approximation for window function computation and 
illustrate this method by computing the Python V window function. 12 We find that the flat-sky 
approximation works well at zero-lag and for pixels at small constant-elevation separation. At 
larger separations the Wigner method is more accurate than the flat-sky approximation. 

The methods developed in this paper are easily extended to other cases not explicitly considered 
here (such as a non-circular non-Gaussian beam, a beam that rotates on the sky, etc.). Elsewhere 
we summarize an analysis of the Python V data that makes use of these methods. 

We acknowledge very valuable discussions with K. Coble, advice from S. Dodelson, M. Drago- 
van, K. Ganga, L. Knox, and the rest of the Python collaboration, and support from NSF CAREER 
grant AST-9875031. 



A. Elliptical Gaussian Beam Function: Normalization & Transforms 

Python V is an example of an experiment with a compact elliptical Gaussian beam function 
(Coble et al. 1999). In such a case it is possible to obtain accurate and useful semi-analytic 
expressions for the Fourier transform B(k) and the spherical harmonic transform bg m (*f-p) of the 
beam function. 

An elliptical Gaussian beam function that is compact enough can be expressed, in a locally 
flat-sky coordinate system (around the beam pointing direction), as 



X\ 2 X2 2 1 



(Al) 



2a{ 2 2a 2 2 

Here x = (x\, x 2 ) are locally flat-sky cartesian coordinates and o\ and a 2 are the beamwidths in 



2 Our flat-sky approximation differs from that implemented in Coble et al. (1999) and Coble (1999). 
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the xi and X2 directions. 13 It is straightforward to establish that this is normalized so that 

/oo />oo 
dxi / dx 2 B(x) = 1 . (A2) 
-oo J — oo 

We shall have need for the expression for the beam in local spherical polar coordinates 7, = 
(6, 4>), around the pointing direction which, without loss of generality, can be assumed to be the 
North Pole 7 P . To fix the orientation of the beam on the sky we choose x\ to lie along = 0. 
Using the local mapping, x\ = Ocoscj) and xi = 6>sin0, we can write 

Q2 



#(7p,7) = 



-exp 



27rai(7 2 ^ [. 2(j2 (^) 
where the "beamwidth" is a function of the polar angle 

2 



a\<t>) 



and a non-circularity parameter 



1 + e 2 shr 



2 2 

V\ - (72 



(A3) 



(A4) 



(A5) 



With the usual accurate approximation, allowed by the rapid fall off of the Gaussian in 9 in eq. (A3), 
it is straightforward to establish that the normalization condition 



/■27T pir 

/ d<p d0sin0B(7p,7) 
Jo Jo 



(A6) 



is satisfied for a(4>) <C 1. 



In the flat-sky approximation the elementary window function defined in eq. (19) depends on 
the Fourier transform of the beam function -B(k). For the elliptical Gaussian beam function of 
eq. (Al) we find 

&iV 2 fc 2 W" 



B(k) = exp 



(A7) 



2 2 

One great advantage of the flat-sky approximation is that Fast Fourier Transform techniques may 
be used to rapidly compute -B(k) for any beam function. 

To use eq. (31) to compute the curved-sky window function we need to compute the spherical 
harmonic transform &£ m (7 P ) of the beam function pointed at the North Pole 7 P . An elliptical 
Gaussian beam function results in a semi-analytic expression that requires numerical evaluation of 
only a single integral. A more complex beam function could require a complete numerical analysis. 



We first note that it is straightforward to show that 

[«7)r = (-irv m (7) 



(A8) 



3 The beam function is also a function of the pointing direction xq, B(x , x), taken to be at the origin here. 
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From eqs. (A3) and (12) we find 



W7p) = T^\ r* n u + m l I He sine e-e 2 /^Pr (cos9) , (A9) 

27T(TiCJ2 y 4vr (e-m)\ J J 

where we have used the usual expression for Y e * m in terms of the associated Legendre function P^ m . 
For a compact beam, and for large £, the 6> integral may be performed using the usual small-0 
approximation , 

lim^ oo £ m p- m (cos0) = J m ([£ + O.5]0) , (A10) 

e.g., eq. (8.722.2) of Gradshteyn & Ryzhik (1994). Using eq. (6.631.7) of Gradshteyn & Ryzhik 
(1994) we find 



, , ^ i- £ y m Ki + m)\ 1 , A11 , 

ir(l + O^Y' 2 y {£ — my. a\a 2 

x / d<t>e- im +f(<i>)e-f W [J (m _ 1)/2 (/ 2 (0)) - / (m+1)/2 (/ 2 (0))] , 
■/ o 

where is the modified Bessel function and 

/W = (£±0|W«. (A12) 

Equation (All) is valid for m > 0; for m < we use this and eq. (A8) for the 6^ m 's. It is 
straightforward to show that for a circular beam eq. (All) reduces to the well known expression. 

Using the reality condition on the beam, eq. ( A8) , one may show that the d(f> integral in 
eq. (All) may be replaced by a [1 + (— l) _m ] J d(f) integral. Clearly, the 6^ m (x3)'s vanish for odd 
m. 14 Since f((f>) is a function of sin 2 (</>), it is straightforward though tedious to show that the 
imaginary part of e~ im< ^ in eq. (All) leads to an expression that vanishes. Thus we have 



[1 + (-1)-"] (I)-" (t + m)\ 1 
bemilp) = ,{£ + 0.5)3/2 V (73^)7 — 

x r#cos(m^/ 3 (0) e -' 2 W [J (ro _i )/2 (/ 2 (0)) - / (m+ i)/ 2 (/ 2 (^))] . 
•/ o 

For an elliptical Gaussian beam function, this approximate semi-analytic spherical harmonic trans- 
form agrees well with the exact fully numerical transform. 



14 This is true provided the integral does not diverge. In fact it is straightforward to establish that it is integrable. 
We use e lm<i> — cos(m0) + ism(m<j)) and consider the real and imaginary parts separately. Both of these are continuous 
functions of (j> over < <f> < tv and are thus Riemann integrable over this interval (see pp. 42 and 63 of Whittaker & 
Watson 1969). 
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FIGURE CAPTIONS 

Fig. 1. — Contour plots of the zeroth order (eq. [22], left panel) and first order (eq. [23], right 
panel) terms in the non-circularity perturbation expansion of the flat-sky approximation elementary 
window function for an elliptical Gaussian beam function experiment. These are computed for the 
nominal FWHM beamwidths of the Python V experiment, 1?02 in elevation and 0?91 in azimuth. 
They are plotted as a function of dimensionless variables fcx, and the two panels are centered on 
the center of the zeroth order circular beam function window function. As expected, the flat-sky 
window function for the circular beam in the left panel is circularly symmetric. For a fixed value 
of k the first order correction in the right panel must be multiplied by — k 2 (<J2 — of) before being 
added to the zeroth order term. For Python V the higher order terms are small compared to the 
first order term and visually have roughly similar structure. 

Fig. 2. — Coefficients \b£ m (-x' 3 ) bg m i(-x' 3 )\ of the six lowest order terms in the perturbation expansion 
of the Wigner method elementary window function (see eqs. [33] and [31]). These are computed 
for the elliptical Gaussian Python V beam function and are plotted as a function of ^orms- For 
the Python V beam function the ellipticity parameter e = 0.26. Non-circularity corrections are 
important for £orms > 1- Note that the peak shifts to higher values of £ for higher order terms, 
relative to the peak position for lower order coefficients. The shape of the curves are independent 
of crms but depend sensitively on the beam function ellipticity. 

Fig. 3. — Contour plots in the azimuth-declination plane (with azimuth along the horizontal axis) 
of terms in the Wigner method perturbation expansion of the elementary window function for the 
Python V experiment. These are plotted for multipole £ = 100 such that ^ctrms ~ 1- The left panel 
shows the zeroth order isotropic term <ioo(7) an d the right panel shows the first order correction 
term {cos(2a) + cos(2/?)} d^il) (see eq. [33]). Here 7 is the angular separation between the central 
pixel and the pixel at the given azimuthal and declination sky coordinates. 

Fig. 4. — Two sets of the eight equal-modulation exact Python V zero-lag complete window func- 
tions Wf m ' m \-f , 7) . The dark solid and light dotted curves correspond to window functions at the 
two extreme declinations 5 = —63° and — 45?4, respectively. Higher modulation window functions 
peak at progressively larger values of £ and with smaller amplitude. The dashed curve is the exact 
Python V zero-lag elementary window function at 5 = —63°. Note that higher modulation com- 
plete window functions peak at £ a few times larger than the inverse beamwidth (<r^ MS ps 140) 
which is where the non-circularity corrections start to become important (see Fig. 2). For the 
same difference in azimuth the angular separation between two equal declination points ( Acf) cos 5) 
is smaller at larger \5\. Consequently, for the same modulation, the window function at S = —63° 
peaks at a larger multipole than the window function at 8 = — 45?4. 

Fig. 5. — Comparisons between Python V flat bandpower equal-modulation covariance matrix 
elements computed using different approximations. Plotted are the relative difference between 
^ m ' m ^(7i,7j) (eq. [41]) computed using an approximate window function and computed using the 
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exact window function, and normalized by dividing by the flat bandpower C£- ' . The approxi- 
mations considered are the circular beam approximation, ("wigOO" in blue), the three successive 
leading improvements to this in the Wigner method perturbation expansion ("wig20" in green, 
"wig22" in red, and "wig40" in black, see eq. [33]; here each successive improvement includes all 
lower order terms) , and the flat-sky approximation ( "flat" in cyan) . Two curves are shown for each 
approximation, those corresponding to modulation 2 (dashed) and 8 (solid). At large £ the modu- 
lation 8 curves converge to a lower accuracy than the corresponding modulation 2 curves because 
the non-circularity correction is more significant for higher modulations. Upper panels show dif- 
ferences for zero-lag covariance matrix elements at the two extreme declinations: a) 5 = —63° and 
b) 5 = — 45?4. Lower panels show differences for non-zero-lag covariance matrix elements. Panel c) 
corresponds to pixels separated in azimuth by 20° at declination S = —63° and panel d) corresponds 
to two neighboring, equal azimuth, pixels at declinations 5 = —63° and —62° (in panel d the solid 
red curve covers the solid black curve for t < 300). The flat-sky approximation fares well in all 
cases except for pixels separated in declination, where it fails even for small-separation pixel pairs, 
see panel d). Panel d) also shows the enormous improvement over the circular beam function ap- 
proximation (wigOO) achieved by accounting for even just the first order Wigner method correction 
(wig20) . Panel c) highlights the need to compute to a large enough value of I to achieve sufficient 
accuracy, e.g., truncating the circular beam function approximation (wigOO) at intermediate i leads 
to an inaccurate result. 
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